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ABSTRACT 



A statement of a method [1] for obtaining an approx- 
imate lower 100(l-a)% confidence limit on System Reliabil- 
ity is made. A discussion of evaluating the accuracy of 
this method on a digital computer is presented. Following 
this the development of a continuity correction factor is 
developed and the accuracies of this refinement for 
several hypothetical syst esented and discussed. 
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CHAPTER I 



PURPOSE OF THE STUDY 

In past years our advancing technology has created 
systems which are becoming more and more complex. In the 
military this has been particularly true, and thus, has 
produced a great need for a method of predicting the over- 
all reliability of these systems. Many problems in the 
area of reliability have not been solved. One of the 
unsolved problems is the formulation of a method for compu- 
ting a confidence interval on system reliability where the 
reliability of the system is computed as the product of the 
reliabilities of each of the components, and the sizes of 
the test samples of the components are not all equal. If 
the test samples are equal, there exist well-known methods 
for computing a confidence interval. Because of the size 
and complexity of present day systems it is necessary for 
major contractors to have many sub-contractors. They may 
be scattered througout the country and may produce a limited 
supply of each component, of which very few may be used for 
reliability testing. These factors plus the general practice 
of testing before the system is assembled points out that 
having the sample sizes all equal is a rare event. 

An ad hoc procedure for computing an approximate con- 
fidence interval has been developed, based upon Bernoulli 
test data from unequal samples of the various components of 
the system [1] . This technique will be called Method I 
throughout this document. 
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An underlying assumption to this procedure is that the 
distribution of the unbiased estimator of the negative 
natural logarithm of system reliability may be approximated 
by a two-parameter Gamma probability distribution. Conse- 
quently, a discrete set of values are approximated by a 
continuous function. Realizing the inherent error in this 
approximation this Thesis has as its purpose 

(1) to develop a correction factor to reduce this 
error; and 

(2) to simulate the distribution of the lower con- 
fidence limit of various systems and examine 
the results achieved with and without this 
correction . 

The fruitfulness of the investigation that follows is 

demonstrated by a discussion of the results. In later 

chapters, it will be explained that the a th percentile 

of the simulated distribution of the lower confidence limit, 

called R _ , . , and the true reliability of the system, R , 

S / i-i V Ot / s 

should be the same; therefore, the difference will be used 
as a measure of the accuracy of Method I. Consider the 
following example: a system is made up of 15 components; 

14 have a sample size of 20 to be tested; and each a failure 
probability of .995; one has a sample size of 150 and prob- 
ability of failure .850. The true reliability, R , is.. 7924. 

5 

The a th percentile point of the distribution of the lower 
100(l-a)% confidence limit, without correction were computed 
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to be: 



a 


R s ,L (a) 


| R s,L(a) R s | 


.05 


.8170 


.0246 


.10 


.8175 


.0251 


.20 


.8000 


.0076 


by adding 


the correction factor 


the results became 


a 


s ,L (a) 


1 R T , ,-R | 
1 s.L(a) si 


.05 


.8113 


.0190 


.10 


.8135 


.0211 


.20 


.7969 


.0045 


The above 


example certainly did 


not have the best 



accuracy of the various systems chosen nor was it the 
worst. The following summary shows the average and vari- 
ance of the error in all the cases studied. 

Without correction: 



a 


Average Error 


Variance of Error 


05 


.0438 


.00213 


10 


.0305 


.00076 


20 


.0109 


.00016 



With correction: 



a 


Average Error 


Variance of Error 


05 


.0286 


.00041 


10 


.0168 


.00015 


20 


.0056 


.00005 



The next chapter plus the appendices contains a com- 
plete explanation of Method I and discusses some further 
approximations made in the method. 

Chapter three discusses the digital computer's role in 
obtaining the distribution of the lower confidence limit of 
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the reliability by Method I. Also mentioned is the measure 
of accuracy used to evaluate the method's precision. 

Chapter four outlines the reason for the correction 
factor and discusses the logic behind its development and 
application to this method. 

In Chapter five the results of various computations 
are presented and discussed. 

The last chapter summarizes the study and states the 
conclusions . 
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CHAPTER II 



EXPLANATION OF METHOD I [1] 

If a system consists of k components in logical series, 
the true system reliability, R g , may be expressed as follows 



R = n 
S i=l 



(1) 



. th 



where p^^ is the true reliability of the i — component. The 

component may be either continuous-operating or cycle type. 

It is desired to find a lower 100(l-a)% confidence 

interval for R . That is, a statistic R T , \ such that 
s ' s,L(a) 

whatever the actual values of p^, P 2 »....P] C , 

PtR s 1 R s ,L (a) ] = 1_a 

The method proposed is to put n^ items of component i 
on test, i=l, 2, ...., k; under the environmental conditions 

defined in the mission and let each operate until failure 
or the mission time is reached, whichever occurs first. 
Letting 

f^ = number of components of type i 



We define: 



and 



that didn't complete their mission, 
p. = 1-q. where q . = ^ i 



(2) 



n . 

l 



k 

= n 

i=l 



The random variable =lnR is now used to obtain the 

s 

lower confidence bound on R g . This is accomplished by 
approximating the distribution of -InR by a two parameter 
gamma distribution and obtaining the confidence bound via 
the gamma distribution. 
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To continue the development define: 



S = -InR = - l In (1-q. ) 

s i=l 1 



( 3 ) 



and, expanding the natural logarithm in an infinite series 



k r 

s = - l 

i=l 



/ \ 1 / \ 2 1 , ,3 

(-q i ) "2 (-q i ) "3 ( " q i ) 



= l l 

i=l j=l 



00 -J 

q J 

i 



(4) 



If each is small, say, .15 or less, the above series 
can be approximated by the first two terms of the series, 
called for ease of expression. That is, 




In Appendix III [ 4 J , it is shown that the error due to the 
aoove truncation is quite small; in fact; 



k 

1 l 

i=l 



T, + l 



iii 3 ^ir 



(6) 



In Appendix I [4] , an unbiased estimator for T^ is devel 
oped and is called T^ . 



T. = 



(7) 



where 

2n^-3 

a i = 2 (n ± -l) 



( 8 ) 



b. 



l 



n . 

l 



nT-1 

l 



( 9 ) 



and as in (2) 
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An unbiased estimator is considered important here because 



the T. are to be added to obtain: 
1 



A p A 

s = y t. 
i=i 1 



(10) 



A 

Therefore, a positive or negative bias on all of the 
would be undesirable. It also follows that S will be an 
unbiased estimator of S. 

An approximate value for the variance of S can be 
shown to be (see Appendices II and IV of [4]) 



k k T . 

Var (S) = l Var (T ± ) = J 



(ID 



i=l J ' i=l -i 

A 

The distribution of S is now fitted with a two param- 
eter gamma family. That is, it is assumed the probability 
distribution of S is given by the density function fg 

(x,r,0) defined as: 

1 



f s (x,r,0) = 



^ x 

X r -Vi3 



, x>0 , r>0 , 0>O 



r U)Q 

0 elsewhere. 



Estimates of r and 0 are found by the Method of Moments 

/\ 

2S 2 

and the fact that is a x„ random variable is used to 

0 A 2r 

obtain a confidence interval for -InR and, thus, from this 

s 

the confidence interval for R can be derived. 

s 



Since S is assumed to be gamma distributed: 
E[S] = 0r 
Var [S] = 0 2 r 
But also from (10) that 



(12) 

( 13 ) 



E[S] = l 



i=l 



T. 

l 
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and 



k T . 
Var[S] = l -i 



i n • 

i=l 1 



Thus, solving the above equations for r and 0, the follow- 
ing expressions are obtained: 



r = 



0 = 




(14) 



(15) 



Thus, an estimator £, for r is taken as 





(16) 



Now to obtain the confidence interval the following 

A 

2 

2r 



2s 

procedure is used. Since — ^ is a Xp- random variable and 



0r = -lnR s> the following probability statements can be 



made : 



1-a = P 



'2S > 2 1 

0 — x l-a,2rJ 



0r < 2S • 



x l-a , 2r 



-InR < 2S • -=r 
S 2 



-InR < 2S 
s — 



x l-ct,2r 

r 



(17) 



x l-a , 2r 

It has been noticed that this last statement has an 

important property; namely 

2r 

2 

x l-a,2r (18) 

is almost independent of r. As an example if a = .10 for 
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2 2 
r = 6 the ratio is ^ and for r = 12 the ratio is 3 — g* 

It is observed that r can be varied by a factor of 2 and 
still have nearly the same value. This effect is even 
less for larger values of r. This suggests that the prob- 
ability distribution of 



[2r] (19) where [2r] denotes- the smallest 

2 

Xl-a [ 2 r] integer greater than or equal to 2 r 

has a very small variance. Therefore, (18) could be re- 
placed by (19) in (17) and, thus, have the following state- 
ment 



1- a = P 



-InR < 
s — 



S 



X 



[2r ] 

2 

1-a , [ 2r ] _ 



( 20 ) 



( 20 ) can be rewritten as 

~-S [ 2 r] 



1-a = P } R g >_ exp 



- X i- a/ [ 2 r] - 



1 

■ 



( 21 ) 



That is a 100 (1-a) % lower confidence limit for R is 

s 



R T / \ 
s ,L (a) 



exp 



~ -S[2f] “ 

2 

- x l-a , [ 2 r ] - 



( 22 ) 
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CHAPTER III 



SIMULATED DISTRIBUTION OF THE CONFIDENCE LIMIT 
In order to evaluate the accuracy of the procedure, 
a digital computer is used to simulate the distribution 
of the S's. The following parameters are used as inputs 
to the computer program: k; n^, i = 1 ,...., k; and p^ = 1-q 
i = 1 , . . . . , k 

In order to have an R g to compare with the simulation the 
program first computes an R from 

b 



k 




The program next simulates values. A three digit ran- 
dom number is obtained using a uniform random number gener- 
ator subroutine. If the random number generated is greater 
than p^, a failure is "counted" by the computer; if the 
random number is less than or equal to p^ t no "count" is 
made. This random number generation is done n^ times for 
each p^ . Thus, the number of failures counted divided by 
the number of units of that component becomes an estimator 
for the unreliability of that component. 



f . 
1 




This process is repeated for all q^; thus giving an estima- 
tor for each q^. Once these estimators are simulated the 
following arithmetic operations are performed by the com- 
puter and stored for further use: 



( 10 ) 



S = l 

i=l 



T. 

1 



f = (I T.) 




I 



T. 

l 

n . 

l 



(23) 



In order to complete the simulation of the distribution 

A 

of R T / \ f replications are made of the above procedure. 

This gives 500 values each of (10) and (23) which are used 
to compute 500 values of 

- ln V(c) " gl2 V X W2f) ,24) 

for each value of a. These 500 points are then sorted, by 
size, by a separate subroutine. The computer now picks out 

the a th percentile of this distribution, call it A. Then 

-A 

R T , > = e is computed and printed out for comparison 

with the true system reliability, R . 

s 

To develop a measure of accuracy for the simulation 
consider the underlying meaning of the probability state- 
ment 

PtR s i VlIc ,) 1 - 1 - “ < 25 » 

If the probability density function of R T , * is plotted 

on a coordinate system, and R g is taken to be a point on 

the abscissa, then 100(l-a)% of the area under the curve 

lies above and a% lies below R . Therefore, R should be 

s s 

equal to the a th percentile of the distribution of R _ , > 

which we have called exp (-A) . 

Thus, |R - exp (-A) | (26) is a measure of the 

s 



accuracy of the procedure under investigation. 
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CHAPTER IV 



DEVELOPMENT OF CONTINUITY CORRECTION FACTOR 
In the previous discussion S was assumed to be a gamma 
distributed random variable. Since a digital computer 
simulation is used to generate the values of S it is real- 
ized that the actual distribution is discrete. Therefore, 
the program is estimating a gamma distribution. As in the 
case of estimating a binomial with a normal distribution 
a more accurate result can be had if the discrete values 
are slightly shifted. 

The following is the development of a correction 
factor (commonly called a continuity correction factor) . 
Consider the following: 



k 



s = 


l T. 

i-1 1 






/s 


a . f i + b i ( f : L | 
i n L 2 \n i , 


2 




T. = 

l 






where a. 

i 


= 2n^-3 and 


b. = n. 
l i 


• 




2 (n^-1) 


n. -1 

l 





The smallest difference between two values of S will 
be the smallest change in one of the T^ terms. For a T^ 
the smallest change is when f^ is increased by one for that 
respective n^ that is the max (n^) 

i=l , 2 , . . . .k 

Let j be that i 3 

n . = max (n . ) 

D i 

i = 1 , 2 , . . .k 
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and in T. increase f. by one to obtain T. ; that is 



(f.+l) b./f.+l 



T. = a. 



J. 



-1UL 



n 



2 in. 

3 



Define the correction factor as 
C 1 



where 



T . = 

3 



Tj* = 



— (T * - 
2 U j 


A 

V 




2n .-3 
3 


f . n . 

. _1 + _1 


f 2 

■ +1 
n . 

3 


2 (n.-l) 


n . 2 (n . -1) 

3 3 


(2n.-3) 


f .+f } 
J 3 




2n^ (n 


r 1 ' 




(2n^-3) 


(f,+l) 

. • 3 A. 


n . 

3 



2 (n_. -1) 



n . 

3 



2 (n_.-l) 



(f j+1) 

2 



n . 

3 



( 28 ) 



(29) 



(30) 



f. + (2n.-l) f . +2 (n . -1) 

= 3 3 3 - 3 

2n . (n.-l) 

3 3 



(31) 



Then 



C, = ± 



[f ^ +(2n j -l)f n +2(n j -l) ] 



[(2n.-3)f +f n ] 
j n n 



2n . (n.-l) 
3 3 



2n . + 2f . -2 

= __J 2 

4n . (n.-l) 

3 3 

n.+f -1 

= _J D_ 

2n j 2 (n j -1) 



(32) 



The corrected S is now defined to be 

A A 

s* = s+c. 



(33) 
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A 

It can be seen from the above calculations that S* 
will have a value even if all the T^ ' s turn out to be zero 
Keeping this in mind and looking at the equation for r 
in equation (23) 




it is observed that in this case r becomes infinite when 
the denominator is zero. Obviously this doesn't make 
sense; therefore, a correction factor is added to the de- 
nominator. Because this correction factor has little in- 
tuitive basis other than keeping f from becoming infinite, 
it was picked quite arbitrarily. The following correction 



factor was chosen: 

A A 

T . - T . 

C~ = -J 1 



n . 

3 



(34) 



and added to 7 — giving a new r ie. £* . This factor 

i-1 n i 

proved to be adequate and didn't appear to bias the final 
results; at least, not in the region of accuracy desired 
in the model. 

The two corrected values are now combined in the prob 
ability statement of equation (21) 

l-o = P [ R s > exp (-S* [2r*]/ Xl ! a _ [2f*l|] < 35 > 

giving the associated lower confidence limit 



R s,L(o) - ex P j-S*[2r*]/: Xl ^ 
which was simulated in the computer. 



[2r*] ) 



(36) 
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CHAPTER V 



RESULTS OF SIMULATIONS 



There were thirteen different combinations of param- 
eters used, representing as many hypothetical systems. 

Table I lists the input parameters and the results. Col- 
umn (1) is self-explanatory. Column (2) is the number of 
components in the system. Colums (3) gives the size of 
the sample tested for each respective component. Column 
(4) is the probability of success for each component of 
the system. Column (5) is the system reliability as com- 
puted directly from the input parameters (2) , (3) and (4) . 

Column (6) lists the three values of a for which a lower 
confidence limit was computed. Column (7) is the exponen- 
tial power of the a th percentile of the distribution. 

The first four cases in Table I are examples of systems 
with different component sample sizes. In Cases 5 through 
13, the sample size is constant in each case. These last 
cases show the effect that system size, sample size and 
probabilities have on the accuracy of the method. The 
effects noted will be discussed after Table II is presented. 

A further check was made on the method by computing 
Case 10 four times, each time starting the random number 
generator at a different point. The results of this test 
showed there wasn*t any difference in the first four decimal 
places of the confidence limit estimate. From this, it is 
concluded the 500 replications of each case is enough to 
smooth out any fluctuations in the generated numbers. 
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TABLE I 



RESULTS OF COMPUTER SIMULATION WITHOUT 
CONTINUITY CORRECTION FACTOR 



(1) 

ZASE 


(2) 

k 


(3) 
n . 

l 






(6)' 

a 


(7) 

exp (-A) 






150 


19 


.995 


.900 




.05 


.7185 






90 


5 


.985 


.980 












75 


125 


.979 


.995 




.10 


.7270 


1 


13 


100 


63 


.988 


.970 


.7233 










125 


125 


.982 


.995 




.20 


.7292 






18 


59 


.980 


.968 












28 




.967 














250 


30 








.05 


.8610 






40 


20 


.990 




.8601 


.10 


.8626 






120 


75 


i = l , — 


-,15 




.20 


.8606 


2 


15 


15 


100 
















130 


90 
















65 


60 
















70 


60 
















130 


















20 


150 


.995 






.05 


.7270 






20 


20 


i-1 , — 


“ , 14 


.7920 


.10 


.7480 






20 


20 








.20 


.7891 


3 


15 


20 


20 


.850 














150 


20 


i=15 














150 


20 
















150 


20 
















150 


















20 i=l 


,-,14 


.995 






.05 


. 8170 


4 


15 






i=l, — 


- , 14 


.7924 


.10 


.8175 






150 


i=15 








.20 


.8000 










.850 


















i=15 
























.05 


1.0000 


5 


5 


50 




.99 




.9510 


.10 


.8958 










i=l , — 


-, 5 




.20 


.9142 
















.05 


.8359 


6 


10 


50 




.99 




.9044 


.10 


.8604 










i=l , — 


- , 10 




.20 


. 8894 
















.05 


.7985 


7 


15 


50 




.99 




.8601 


.10 


.8325 










1=1 , — 


~ , 15 
! 




.20 : 


.8506 
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TABLE I (Continued) 



RESULTS OF COMPUTER SIMULATION WITHOUT 
CONTINUITY CORRECTION FACTOR 



(1) 

CASE 


(2) 

k 


(3) ” 
n . 

i 


— _ 

Pi 




Te y ■" 

a 


T7) 

exp (-A) 


8 


30 


20 


.99 

i=l , ,30 


.7397 


.05 

.10 

.20 


.5767 

.6647 

.7059 


9 


5 


100 


.95 

i=l , r 5 


.7738 


.05 

.10 

.20 


.7636 

.7652 

.7696 


10 


15 


20 


.995 

i=l, ,14 

.85 

i=15 


.7924 


.05 

.10 

.20 


.6862 

.7048 

.7666 


11 


15 


40 


.995 

i=l, ,14 

.85 

i=15 


.7924 


.05 

.10 

.20 


.7683 

.7825 

.7935 


12 


15 


50 


.995 

i=l, ,14 

.85 

i=15 


.7924 


.05 

.10 

.20 


.7704 

.7829 

.7925 


13 

l ■ 


15 


100 


.995 

i=l, — ,14 

.85 

i=15 


.7924 


.05 

.10 

.20 

j 


.7780 

.7885 

.7918 
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Table II lists the results of the cases with and with- 



out the correction factor added. Columns (1), (2) and (3) 

were described before. Column (4) is the exponential 
power of the a th percentile of the distribution with the 
correction added. Column (5) shows the value of the meas- 
ure of accuracy used in evaluating Method I with correction 
added. Column (6) is the same as (5) but without the 
correction added. 

From Table II, it is quite apparent, on the whole, 
that the error was reduced by the addition of the correction 
factor. It also appears the correction has less effect on 
cases where the sample sizes (n^) vary [Cases: 1, 2, 3, 4], 
but still, obviously, worth the effort. Cases 5, 6 and 7 
seemed to indicate the number of components (k) in the 
system has no effect on accuracy. Cases 10, 11, 12 and 13 
appear to point out that as the sample size (n.) increases 

\V\C < 

the accuracy versus probabilities. This concludes the 
discussion of the most observable facts that can be viewed 
from the tables. 




3 4tet> GL&AiiJl Wj 
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TABLE II 



COMPARISON OF METHOD I WITH (*) AND 
WITHOUT CONTINUITY CORRECTION FACTOR 





(1) 


(2) 


(T> 


(4) 


(5) 




.(6) 


3ASE 


R s 


a 


exp (-A*) 


exp ( —A* ) -R s 




exp (-A) -R s 






.05 


.7184 


.0049 




.0048 


1 


.7233 


.10 


.7270 


.0037 




.0037 






.20 


.7291 


.0058 




.0038 






.05 


.8610 


.0010 




.0010 


2 


.8601 


.10 


.8626 


.0025 




.0026 






.20 


.8606 


.0006 




.0006 






.05 


.7219 


.0705 




.0654 


3 


.7924 


.10 


.7459 


.0465 




.0444 






.20 


.7873 


.0051 




.0033 






.05 


.8113 


.0190 




.0246 


4 


. 7924 


.10 


.8135 


.0211 




.0251 






.20 


.7969 


.0045 




.0076 






.05 


. 9961 


.0039 




.0490 


5 


.9510 


.10 


.9269 


.0241 




.0552 






.20 


.9522 


.0012 




.0368 






.05 


. 8628 


.0415 


1 


.0685 


6 


.9044 


i .10 


.8953 


.0090 


! 

L 


.0440 




i 


.20 


.8998 


.0045 


f 

L. 


.0150 






.05 


.8384 


.0216 


! 


.0616 


7 


.8601 


.10 


.8480 


.0120 


i 


.0276 






.20 


.8567 


.0033 


H 

; 


.0095 






.05 


.6830 


.0567 


i 

I 

I 


.1630 


8 


.7397 


.10 


.7070 


.0327 




.0750 






.20 


.7220 


.0177 




.0338 




1 

\ 


.05 


.7635 


.0103 




.0102 


9 


| .7738 


.10 


.7651 


.0086 


i 


.0085 




» 


.20 


.7696 


.0042 




.0042 






.05 


.7495 


.0429 




.1062 


10 


.7924 


.10 


.7591 


.0330 


t 

i 


.0876 




| 

\ 


.20 


.7684 


.0240 

! ! 

i ; 


i 


.0258 
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TABLE II (Continued) 



COMPARISON OF METHOD I WITH (*) AND 
WITHOUT CONTINUITY CORRECTION FACTOR 





^1) 


[21 


(3) 


(4) 


( 5 ) 


(61 


2ASE 


R s 


a 


exp (-A* ) 


exp (-A*)-R s 


exp (A-A) -R s 






.05 


.7714 


.0210 


.0241 


11 


.7924 


.10 


.7820 


.0104 


.0099 






.20 


.7931 


.0007 


.0011 






.05 


.7700 


.0224 


.0220 


12 


.7924 


.10 


.7826 


.0098 


.0095 






.20 


.7923 


.0001 


.0001 






.05 


.7779 


.0145 


.0144 


13 


.7924 


.10 


.7884 


.0040 


.0039 






.20 


.7917 


.0007 


.0006 
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CHAPTER VI 



SUMMARY AND CONCLUSIONS 



In the preceding four chapters a procedure [1] has 
been explained which will approximate a lower 100(l-a)% 
confidence interval on system reliability. The simulation 
of this method on a computer has been explained and the 
technique used to measure the accuracy of the method was 
stated. A continuity correction factor has been developed 
which enables the discrete computations of a digital com- 
puter to more closely approximate a continuous distribution. 

The simulation was run for thirteen cases with and 
without the continuity correction factor. These were not 
an adequate amount of cases on which to base any concrete 
conclusions. It is felt these cases are sufficient to 
show the merit of the correction factor. The following 
summary as stated in Chapter I demonstrates the validity 
of the procedure. 

Without Correction: 

a Average Error Variance of Error 



05 


.0438 


10 


.0305 


20 


.0109 



.00213 
.00076 
. 00016 



With Correction: 



a 



Average Error 



Variance of Error 



05 


.0286 


10 


.0168 


20 


.0056 



.00041 

.00015 

.00005 



From Tables I and II in Chapter V, it can be seen 
the average error for those cases where sample size (n^) 
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is varied would be less than shown above. It must then 



be concluded that Method I with the correction factor 
will be an effective method for approximating the lower 
100(1 - a) % confidence limit for system reliability. 

It should be pointed out that another study; ob- 
viously, could be done in this same area. By varying the 
input parameters (k, n^ , p^) and determining their 
effects, one could possibly establish limits on their 
variation to keep within a desired degree of accuracy. 

In closing, I would like to express my appreciation 
to Dr. W. Max Woods of the Operations Analysis Department 
for pointing out this problem and serving as a very 
willing advisor. I also wish to extend my thanks to Dr. 
Rex H. Shudde for supplying the basic computer program 
for Method I which required slight alteration to meet 
my needs . 
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APPENDIX I 



DERIVATION OF T . , THE UNBIASED ESTIMATOR OF T. 

l ' l 



In order that T. be an unbiased estimator of T, , 

1 l 



it is necessary that 



T. = E [T. ] 

l l 



Rewriting the above expression, using (5) and (7) 



*i + -T- = E 



a . f . b . f . 

1 1 + 1 i 



2 -i 



n . 

l 



and, substituting (2) 



*i + — 



= E 



a . f . b.f. 

l l li 



2 -i 



n . . 2 

l 2n . 

l 



a . b . ~ 

— — ■ E(f.) + — E[f, 1 

n i 1 2n. 2 1 

1 



a . n . q . b . 



n . 0 2 

l 2n. 

l 



, z z z , . 

(n. q. -n.q. +n.q.) 

l M i i^i i^i 



b-q. , 

a . q . + 1 ~ 1 ( 1 - — ) + b . q . 

i M i 2 n . i^i 

l 



b. b. q.2 

(as + - + (b . 

l ^ l 
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Then equating coefficients of 




a 



i 



+ 



b. 

1 __ 

2n . 

1 



1 




and solving simultaneously 








X ( b i) 










Substituting back 



=> 




b i) 




1 



3 



b. + — 
1 n i 



2 



b. = 

l 



n . 



i 



n . - 



i 



1 
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and again substituting 



= — (3 
2 



n . 



- ( 



n , 



- 1 



) ) 



2n . - 3 

l 

2 (n i - 1) 



Thus , 



l 



2n . - 3 n . 

2 (n^ - 1) ^i + n^ - 1 




is an unbiased estimator of T.. 

l 
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APPENDIX II 



COMPUTATION OF VARIANCE OF S 



I k ^ \ k 

Var (S) = Var l T . = £ Var (T.) 

\i=l 1 I i=l 1 



Var (T i ) = E[T i 2 ] - E 2 [T ± ] 



= E [T. 2 ] - T . 2 
1 1 



E[T i 2 ] = E[(a.q. + - ) 2 ] 



,2. 4 

_ pr 2^2, ^ /s 3 , i n 

= E [ a . q . + a . b . q . + 7 ] 

1^1 1 1^1 4 



2 u K 2 
a . 9 a . b . o b . . 

= -iy E[f. 2 ] + -±g± E[f. 3 ] + E[f. 4 ] 



n . 
1 



n . 
1 



4n i 



From [ 3 ] 



2 , 2 2 2 

E [ f . ] = n . q . - n . q . +n.q. 

1 1 1^1 1^1 



E[f i J ] = E [ (f i - E[f i ]) 3 ] + 3E[f i ]E[f i 2 ] - 2E J [f i ] 



33 , 2 3 0 3 ,22 , 2 ^ 

= n. q. - 3n. q. + 2n.q. + 3nrq. - 3n.q. + n.q. 

1^1 1^1 1^1 1^1 1^1 li 



E[f i 4 ] = E[(£. - E [f i ] ) 4 ] + 4E[f.]E[f i 3 ] - 6E 2 [f ± ] E [f ± 2 ] + 



3E [f.] 
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44 , 3 3 ^ 2 4 , „ 2 3 ^ 2 2 

n i <Ji - 6n i <Ji + lln i "Ji - 18n i ■>! + qi 



. 4 , - 3 2 

6n.q. + 12n.q. - 7n.q. + n.q. 

ii i^i ii li 



Therefore 



- (2n. - 3) oo o 

E[ *i ] = 7 2 ~7T ~ 2 2~ {n i - n A + n i*i } + 

z ( n . - 1 ) n . 

l i 



(2n^ -3 ) n i33 23 .23 

— 772— J {n i ‘•i ' 3n i ‘Ji + 3n i ‘Si + 

2 (n . - 1) n . 

1 1 



~ 2 2 ^ 2 

1 ^1 1^1 1^1 



n i > r 4 4 3 4,>- 3 3,,, 2 4 

in. q- “ 6n. q. + 6n . q. + lln. q. 

, v 2 4 1 ^1 1 ^1 1 1 ^1 



4 (n . - 1) n 
1 1 



to 2 3 , - 2 2 . 4 . 0 3 - 2 

18n. q. + 7n. q. - 6n.q. + 12n.q. - 7n.q. + 

1 1 i 4 i 1^1 1^1 



n . q . ^ 
1 ^ 1 ) 



^ 2 

E[T i ^] 



(4n. 2 - 12n . + 9) (n.q. 2 - q . 2 + q . ) + (4n. - 6) 

JL -L ^ X. «L 

4 (n . - 1) 2 n . 

1 1 



,23 r> 3 , „ 3 , 0 2 _ 2 J , 

(n i q ± - 3n iq± + 2 qi + 3n iqi - 3 qi + qi ) + 

4 (n^ - 1) 2 n^ 
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(n . q . - 6n . q . + 

v l l 



6n . 



+ lln.q. - 18n.q. + 



l 



l 



4 (n . - 1) n . 

l i 



7 n . qi 2 - 6q^ 4 + 1^*3^ 3 - 7q. 2 + q. ) 

4(n. - 1) 2 n . 

1 1 

3 4 3 2 2 4 3 

n. (q. + 4q. + 4q. ) + n- (-6q, - 12q. 

_ i 2i l l 

4n i (n i - l) 2 

4q i 2 + 4q ± ) + n^llq^ + 8q i 3 - 2 qi 2 - 8q i ) + 

4n i (n i - l) 2 



(-6q i 4 + 2q ± 2 + 4q A ) 
4n i (n i - l) 2 



Let 



4 3 2 

A = q i + 4q ± + 4 qi 



B = -6q ± 4 - 12q i 3 - 4 qi 2 + 4 qi 



C = ll qi 4 + 8q. 3 - 2q. 2 - 8q.. 



D = -6q i 4 + 2q ± 2 + 4q ± 
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Then 



„ An . 3 + Bn . 2 + Cn . + D 
E[t 2 j = — ±_j — i — 

r 4n. - 8n . + 4n . 

111 



A , B + 2A , <C + 2B + 3A)n i + (D - B - 2A » 

T + " Tn" + J ~T2 

i 4n. (n. - 1) 

l i 



4 3 2 /0 4 . 3 0 2v 

2 ~ q i ~ q i + _Ji + q i (2q i - 4q. + 2q. ) n ■ 

1 n i 4n.(n. - l) 2 



-(2q i 4 - 4 qi 3 + 2 qi 2 ) 

4 n i ( n i " 1)2 



2 _ 2 ,_ 3 . _ 4 W 2,_ 1} 2 



2 q A q ± q A (q A + q i >* q i ( q i 

= T. + — + i ' + - . ---r-i 

l n . 2n . 2n . n . 2n . (n . 

ill ‘ 



2n i ( n i " 1) 



a O 

E[T.^] 



= VarfT^ 



T. q . 

-i + -±- 

n i n i 



1 _ (q, + qi ) + 

2 1 1 



(q A - i) 

2(n i - 1) 



which, when simplified, gives 



T . q 

Var ($. ) = — + — 



n i n i 



3 q . 2 - n i (2 q . 2 + 2 q . - 1) 
_ 1} 



t ; 

1 + 



> n d r 2(n 4 - 1) 



q j 
t . 

i 



3q, 2 



1 + — i 2q. - 2 q . 

n. ^i ^i 
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which we approximate by 



T. 

Var(T i ) =‘ji 
1 



Thus 



Var (§) 




T. 

1 



n . 



and the truncation error is 



k 

I 

i=l 



2 (n ± - 1) 



1 + 



3q. 



n . 

i 



- 2q A - 2q ± 



This error is further discussed in Appendix IV. 
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APPENDIX III 



ANALYSIS OF THE ERROR DUE TO TRUNCATION 



In the transitions from (4) to (5) an infinite 
series was truncated to its first two terms. Examine 
the inner summation of (4) 



00 



l 

j=l 





(37) 



= T. + R„ 

l 2 



where R 2 is the remainder of the series after two terms 

th s t 

have been written out. If we test the n 1 " and (n 4- 1) 
terms of this series by the ratio test, 

nq. 
n + 1 



n+ l / L 1 
/ n + 1 

q L / n 



we see that this expression is always less than q, < 1, 
which implies the series is convergent. Then by Theorem 
24, p. 328 of [2] 



lira 



n->-°° 



nq. 

HTT 



which implies that 

R n ± q i n+1 /(n + 1) (1 - r) , r > q. 
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Since r ^ q^, let us choose r = q. ^ 1 so as to make 

R as small as possible 
n c 

R £q. n + 1 / (n+ 1) ( 1 - q . ) 

and evaluating at n = 2 

R 2 * q i / 3(1 " < ^i ) 



Since q. << 1 in most cases of interest, choose q. = 
.15, for example, which is about the largest value it 

_3 

will ever take on, and we see that R 2 i 1.4 x 10 
Substituting equation (38) into (37), 



oo 



l 

j=l 



V 



< T. 




/ 3d - qi) 



and substituting (39) into (4) 



S * I (T + q 3 / 3(1 - q.)) 
i=l 1 1 ' 1 

which upon distributing the summation sign becomes 
equation (6). For q^ = .05, still quite a large value 
for q^ to take on, and k = 15 

k r 

S < y T. + 6 x 10 , 

" i=l 1 



(38) 



(39) 
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so the approximation 



k 

S = l T. 
i=l 

is quite good. 



APPENDIX IV 



INVESTIGATION OF ERROR DUE TO APPROXIMATION IN Var (S) 



In Appendix II the Variance of the random variable S is 



computed, and it is seen that the value used in (11) before 
differs from the true variance by a factor 



Upon embarking upon an investigation of the size of this term 
it is quickly seen that it is a task which, due to the amount 
of work involved, is beyond the scope of this thesis. A few 
general comments on this term are in order, however. 

Being a function of both the and n^ , equation (40) 
may, for certain combinations of these parameters, actually 

A 

be less than zero and reduce the variance of S. In order to 
determine when this occurs, it would be necessary to use a 
digital computer to compute and plot for each set of n^, 
i=l,..., k, a curve over the range of 0 < q. < 1; and, 
similarly, to compute and plot for each set of q^ , i=l,...,k, 
a curve over the range of interest of n., 2 £ n^ £ N, where 
N might be 150. 

If the summand of (40) is plotted as a function of a 
single q^, with the n^ held fixed, the curve is seen to be 
bell-shaped over the range 0 £ q £ 1, but within the real- 
istic range that q^ may be expected to take on 0 £ q^ £ .15, 
(40) is a non-decreasing function, and strictly increasing 
unless a point of inflection occurs. 




2 



2 



(40) 
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If the summand of (40) is plotted as a function of a 

single n^, where n^ is of course interger-valued , it is seen 

to be a monotone-decreasing step-function, asymptotic to the 

n . axis . 

1 

As a crude bound, on the size of the error, consider 
the summand of (40) 



J5i 

2 (n^-l) 



1 + 



n . 



- 2q A - 2qi 



2 

2n . 

l 



1 + 



3q^ 



n . 

i 



- 2qi - 2q £ 



2i 

2n . 

l 



1 + 



n r 
1 



2n . 
i 



1 + 



n . 



which implies that 



2 

*i 


3 9i 2 

I _ 1 _ . _ X _ 


2qi - 2q t 




k 1 

< y 1 




i — 1 

1 

c 

CM 


n . 




2n . 


1 + K7 


1 


i 


J 




i=l l 


1 



which is obviously not a least upper bound, but an example 
is informative. 
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Let n = 20 for all i, and k = 15 



k 



l 

i=l 



1 

2n . 

l 




= .4 



which admittedly is of an undesirable magnitude for an omit- 
ted term, but (40) is undoubtedly much less than .4. 
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